# Analysis script for results reported in 
#Coproduction of Core and Complementary Tasks in Times of Service Decline: Experimental Evidence
#
#
# Setup and packages ####
# Working direction should match the location of the data file
setwd("INSERT HERE") # Working directory

library(tidyverse) # Loading of packages used
library(stargazer)
library(psych)
library(GGally)
library(patchwork)
library(ggsignif)
library(janitor)
library(rstatix)
library(effectsize)


ggplot2::theme_set(ggplot2::theme_minimal()) # Setting the baseline theme in ggplot

# Import data ####

svdat <- read.csv("svdat.csv")

svdat <- svdat %>% 
  mutate(declinecore = as.character(declinecore))

# Manipulations check #####
lm(declbel ~ declcue, svdat) %>% summary()


#Table 1. Descriptive statistics #### 


descvars1 <- svdat %>% dplyr::select(
  cintent,coreintent, compintent
  ) 

desc1 <- describe(descvars1) %>% 
  mutate(mean = round(mean, 3), sd = round(sd, 3))

desc1 <- desc1 %>% add_rownames(var = "var") %>% 
  dplyr::select(
    var, n, mean, sd, min, max
  ) 


descvars2 <- svdat %>% 
  filter(!is.na(compintent)) %>% 
  dplyr::select(
  female, qage, edu, partner, east, 
  ) 
  

desc2 <- describe(descvars2) %>% 
  mutate(mean = round(mean, 3), sd = round(sd, 3))

desc2 <- desc2 %>% add_rownames(var = "var") %>% 
  dplyr::select(
    var, n, mean, sd, min, max
  ) 

descmerge <- bind_rows(desc1, desc2)

descmerge

# Figure 2 ####

copr.mod1 <- lm(coreintent ~ declcue, svdat) 
core.mod1 <- lm(compintent ~ declcue, svdat) 
comp.mod1 <- lm(cintent ~ declcue, svdat)    

fig2ci95 <- bind_rows(   
  tidy(copr.mod1, conf.int = T, conf.level = .95) %>% 
    mutate(dv = "Overall coproduction") %>% 
    mutate(plot ="A") %>% 
    mutate(nobs =nobs(copr.mod1)),
  tidy(core.mod1, conf.int = T, conf.level = .95) %>% 
    mutate(dv = "Core tasks") %>% 
    mutate(plot ="A") %>% 
    mutate(nobs =nobs(core.mod1)),
  tidy(comp.mod1, conf.int = T, conf.level = .95) %>% 
    mutate(dv = "Complementary tasks") %>% 
    mutate(plot ="A") %>%
    mutate(nobs = nobs(comp.mod1))) %>% 
  mutate(level = .95) %>%
  filter(term != "(Intercept)") %>% 
  mutate(axisgroup = case_when(
    dv == "Overall coproduction" ~ 1,
    dv == "Core tasks" ~ 2,
    dv == "Complementary tasks" ~ 3
  )) 

fig2ci90 <- bind_rows(   
  tidy(copr.mod1, conf.int = T, conf.level = .9) %>% 
    mutate(dv = "Overall coproduction") %>% 
    mutate(plot ="A"),
  tidy(core.mod1, conf.int = T, conf.level = .9) %>% 
    mutate(dv = "Core tasks") %>% 
    mutate(plot ="A"),
  tidy(comp.mod1, conf.int = T, conf.level = .9) %>% 
    mutate(dv = "Complementary tasks") %>% 
    mutate(plot ="A")) %>% 
  mutate(level = .9) %>%
  filter(term != "(Intercept)") %>% 
  mutate(axisgroup = case_when(
    dv == "Overall coproduction" ~ 1,
    dv == "Core tasks" ~ 2,
    dv == "Complementary tasks" ~ 3
  )) 

pan_a_title <- expression(paste(bold("Panel A"), "   Unstandardized estimates, OLS"))

fig2panela <- 
  ggplot(data = fig2ci95,  # Plotting for overall and core/comp task
         aes(x=term, y = estimate, group = fct_reorder(as.factor(dv), -axisgroup), ymin= conf.low, ymax = conf.high))+
  coord_flip()+
  geom_hline(yintercept = 0, color = "black", alpha = 0.5, size = 0.5, linetype = "dashed") +
  geom_linerange(position = position_dodge(width=0.9), color="grey40") + # 95% data
  geom_linerange(data = subset(fig2ci90),  #90 % data 
                 aes(term,
                     y = estimate,
                     ymin= conf.low, 
                     ymax = conf.high, 
                     group = fct_reorder(as.factor(dv), -axisgroup)),
                 position= position_dodge(width = 0.9), size = 1.3, color = "grey60")+
  
  geom_point(position = position_dodge(width=0.9))+
  geom_text(aes(y = conf.low, label = dv), 
            position = position_dodge(width = 0.9), hjust= 1.1, 
            color = ifelse(fig2ci95$axisgroup %in% c(1,2,3), "black", "grey40"),
            fontface = ifelse(fig2ci95$axisgroup %in% c(1,2,3), "plain", "plain"),
            size = ifelse(fig2ci95$axisgroup %in% c(1,2,3), 3.2, 3))+
  geom_text(aes(y = estimate, 
                label = gsub("0\\.", "\\.", 
                             paste("p =",
                                   round(p.value, 3)))), 
            position = position_dodge(width = 0.9), vjust= 1.5, 
            size = 3,
            fontface = ifelse(fig2ci95$axisgroup %in% c(1,4,7), "plain", "plain"),
            color = ifelse(fig2ci95$axisgroup %in% c(1,4,7), "grey40", "grey40"))+
  geom_text(aes(y = estimate, 
                label = gsub("0\\.", "\\.", 
                             paste("n =",
                                   round(nobs, 3)))), 
            position = position_dodge(width = 0.9), vjust= 3, 
            size = 3,
            fontface = ifelse(fig2ci95$axisgroup %in% c(1,4,7), "plain", "plain"),
            color = ifelse(fig2ci95$axisgroup %in% c(1,4,7), "grey40", "grey40"))+
  
  
  labs(title = "", 
       x="",
       y=pan_a_title)+
  scale_y_continuous(breaks=c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3),
                     limits = c(-0.3, 0.3))+
  theme_bw()+
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.title.y = element_blank())+
  scale_x_discrete(labels= c(""))

# For panel B 

#First, calculation of cohens d (pooled sd)
cohencopr95 <- cohens_d(svdat$cintent[svdat$declcue == 1], svdat$cintent[svdat$declcue == 0], ci = 0.95)
cohencore95 <- cohens_d(svdat$coreintent[svdat$declcue == 1], svdat$coreintent[svdat$declcue == 0], ci = 0.95)
cohencomp95 <- cohens_d(svdat$compintent[svdat$declcue == 1], svdat$compintent[svdat$declcue == 0], ci = 0.95)

cohencopr90 <- cohens_d(svdat$cintent[svdat$declcue == 1], svdat$cintent[svdat$declcue == 0], ci = 0.90)
cohencore90 <- cohens_d(svdat$coreintent[svdat$declcue == 1], svdat$coreintent[svdat$declcue == 0], ci = 0.90)
cohencomp90 <- cohens_d(svdat$compintent[svdat$declcue == 1], svdat$compintent[svdat$declcue == 0], ci = 0.90)


#Variable ID added
cohencopr95 <- cohencopr95 %>% 
  mutate(dv = "Overall coproduction") %>% 
  data.frame()
cohencopr90 <- cohencopr90 %>%
  mutate(dv = "Overall coproduction") %>% 
  data.frame()

cohencore95 <- cohencore95 %>% 
  mutate(dv = "Core tasks") %>% 
  data.frame()
cohencore90 <- cohencore90 %>% 
  mutate(dv = "Core tasks") %>% 
  data.frame()

cohencomp95 <- cohencomp95 %>% 
  mutate(dv = "Complementary tasks") %>% 
  data.frame()
cohencomp90 <- cohencomp90 %>% 
  mutate(dv = "Complementary tasks") %>% 
  data.frame()

fig2panbci95 <- bind_rows(   
  cohencopr95,
  cohencore95,
  cohencomp95) %>% 
  mutate(axisgroup = case_when(
    dv == "Overall coproduction" ~ 1,
    dv == "Core tasks" ~ 2,
    dv == "Complementary tasks" ~ 3
  )) %>% 
  mutate(term = "declcue")

fig2panbci90 <- bind_rows(   
  cohencopr90,
  cohencore90,
  cohencomp90) %>% 
  mutate(axisgroup = case_when(
    dv == "Overall coproduction" ~ 1,
    dv == "Core tasks" ~ 2,
    dv == "Complementary tasks" ~ 3
  )) %>% 
  mutate(term = "declcue")

pan_b_title <- expression(paste(bold("Panel B"), "   Standardized estimates, Cohen's d"))

fig2panelb <- 
  ggplot(data = fig2panbci95,  # Plotting for overall and core/comp task
         aes(x=term, y = Cohens_d, group = fct_reorder(as.factor(dv), -axisgroup), ymin= CI_low, ymax = CI_high))+
  coord_flip()+
  geom_hline(yintercept = 0, color = "black", alpha = 0.5, size = 0.5, linetype = "dashed") +
  geom_linerange(position = position_dodge(width=0.9), color="grey40") + # 95% data
  geom_linerange(data = subset(fig2panbci90),  #90 % data 
                 aes(term,
                     y = Cohens_d,
                     ymin= CI_low, 
                     ymax = CI_high, 
                     group = fct_reorder(as.factor(dv), -axisgroup)),
                 position= position_dodge(width = 0.9), size = 1.3, color = "grey60")+
  
  geom_point(position = position_dodge(width=0.9))+
  geom_text(aes(y = CI_low, label = dv), 
            position = position_dodge(width = 0.9), hjust= 1.1, 
            color = ifelse(fig2panbci95$axisgroup %in% c(1,2,3), "black", "grey40"),
            fontface = ifelse(fig2panbci95$axisgroup %in% c(1,2,3), "plain", "plain"),
            size = ifelse(fig2panbci95$axisgroup %in% c(1,2,3), 3.2, 3))+
  labs(title = "", 
       x="",
       y=pan_b_title)+
  scale_y_continuous(breaks=c(-0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3),
                     limits = c(-0.3, 0.3))+
  theme_bw()+
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.title.y = element_blank())+
  scale_x_discrete(labels= c(""))

figure2 <- ( fig2panela | fig2panelb )

ggsave(figure2, 
       filename = "figure2.png", dpi = 300, 
       type = "cairo", width = 22, height = 11, units = "cm")

#Figure 3####
  
  # 2.1.1 Panel A
  
#Interaction model

copr.mod2 <- lm(cintent ~ declcue*corecue, svdat) # Coproduction overall
comp.mod2 <- lm(compintent ~ declcue*corecue, svdat) # Core tasks
core.mod2 <- lm(coreintent ~ declcue*corecue, svdat) # Complementary tasks

pan_a_fig3_title <- expression(paste(bold("Panel A"), " Willingness to coproduce core tasks across groups"))
  
  
  source("https://raw.githubusercontent.com/datavizpyr/data/master/half_flat_violinplot.R")
  
  fig3panela <- ggplot(svdat, aes(declinecore,
                    coreintent, 
                    fill=declinecore)) +
    geom_flat_violin(alpha = 0.7) +
    annotate("text", label = "Control \n Complementary Task \n Encouragement", x = 1, y = 5.5, size = 2.9, colour = "black")+
    annotate("text", label = "Control \n Core Task \n Encouragement", x = 2, y = 5.5, size = 2.9, colour = "black")+
    annotate("text", label = "Service Decline Cue \n Complementary Task \n Encouragement", x = 3, y = 5.5, size = 2.9, colour = "black")+
    annotate("text", label = "Service Decline Cue \n Core Task \n Encouragement", x = 4, y = 5.5, size = 2.9, colour = "black")+
    stat_summary(fun.y = mean, color = "white", position = position_dodge(0.75),
                 geom = "point", shape = 21, stroke = 2, size = 3, fill = "black",
                 show.legend = FALSE)+
    theme(legend.position="none") + 
    theme_bw()+
    theme(
      panel.grid.minor = element_blank(),
      panel.grid.major.x = element_blank(),
      axis.text.x = element_blank(),
      axis.ticks.x = element_blank(),
      axis.title.y = element_blank(),
      legend.position = "none")+
    scale_y_continuous(limits = c(1,6), breaks = c(1,2,3,4,5))+
    scale_fill_grey()+
    scale_color_grey()+
    labs(title="", x = pan_a_fig3_title)
  
  
# Panel B 


  fig3ci95 <- bind_rows(   
  tidy(core.mod2, conf.int = T, conf.level = .95)) %>% 
    mutate(nobs = nobs(core.mod2)) %>% 
  mutate(axisgroup = case_when(
    term == "declcue" ~ 1,
    term == "corecue" ~ 2,
    term == "declcue:corecue" ~ 3
  )) %>%  
    mutate(labz = case_when(
      term == "declcue" ~ "Service decline cue",
      term == "corecue" ~ "Core task encouragement",
      term == "declcue:corecue" ~ "Service decline cue × Core task encouragement")) %>%
    filter(term != "(Intercept)")

  fig3ci90 <- bind_rows(   
    tidy(core.mod2, conf.int = T, conf.level = .90)) %>% 
    mutate(nobs = nobs(core.mod2)) %>% 
    mutate(axisgroup = case_when(
      term == "declcue" ~ 1,
      term == "corecue" ~ 2,
      term == "declcue:corecue" ~ 3
    )) %>%  
    mutate(labz = case_when(
      term == "declcue" ~ "Service decline cue",
      term == "corecue" ~ "Core task encouragement",
      term == "declcue:corecue" ~ "Service decline cue × Core task encouragement")) %>% 
    filter(term != "(Intercept)")

pan_b_fig3_title <- expression(paste(bold("Panel B"), " Effects on willingness to coproduce core tasks"))

fig3panelb <- 
  ggplot(data = fig3ci95,  # Plotting for overall and core/comp task
         aes(x=term, y = estimate, group = fct_reorder(as.factor(labz), axisgroup), ymin= conf.low, ymax = conf.high))+
  coord_flip()+
  geom_rect(aes(xmin = 0, 
                ymin = -1*(sd(svdat$coreintent, na.rm=T) * 0.2), 
                xmax = 4, ymax = sd(svdat$coreintent, na.rm=T) * 0.2), fill = "grey80", color = "white", alpha = .07)+
  geom_text(aes(x = 0.3, 
                y = -sd(svdat$coreintent, na.rm =T)*0.2),
            label = "- 0.2 SD",
            size = 3, 
            hjust = .5,
            colour = "grey40"
  )+
  geom_text(aes(x = 0.3, 
                y = sd(svdat$coreintent, na.rm =T)*0.2),
            label = "0.2 SD",
            size = 3, 
            hjust = 0.5,
            colour = "grey40"
  )+
  
  geom_hline(yintercept = 0, color = "black", alpha = 0.5, size = 0.5, linetype = "dashed") +
  geom_linerange(position = position_dodge(width=0.9), color="grey40") + # 95% data
  geom_linerange(data = subset(fig3ci90),  #90 % data 
                 aes(term,
                     y = estimate,
                     ymin= conf.low, 
                     ymax = conf.high, 
                     group = fct_reorder(as.factor(labz), axisgroup)),
                 position= position_dodge(width = 0.9), size = 1.3, color = "grey60")+
  
  geom_point(position = position_dodge(width=0.9))+
  geom_text(aes(y = estimate, label = labz), 
            position = position_dodge(width = 0.9), 
            vjust= -1.5,
            size = 3.2, 
            color = "black")+
  geom_text(aes(y = estimate - 0.065, 
                label = gsub("0\\.", "\\.", 
                             paste("p =",
                                   round(p.value, 3)))), 
            position = position_dodge(width = 0.9), vjust= 1.9, hjust = 0,  
            size = 3,
            fontface = "plain",
            color = "grey40")+
  geom_text(aes(y = estimate - 0.065, 
                label = gsub("0\\.", "\\.", 
                             paste("n =",
                                   round(nobs, 3)))), 
            position = position_dodge(width = 0.9), vjust= 3.4, hjust = 0,
            size = 3,
            fontface = "plain",
            color = "grey40")+
    labs(title = "", 
       x="",
       y=pan_b_fig3_title)+
  scale_y_continuous(breaks=c(-0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6),
                     limits = c(-0.6, 0.6))+
  theme_bw()+
  theme(
    panel.grid.minor = element_blank(),
    panel.grid.major.y = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    axis.title.y = element_blank(),
    legend.position = "none")+
  scale_x_discrete(labels= c(""))

# Combing the panels 

figure3 <- (fig3panela | fig3panelb)

#Exporting figure 3

ggsave(figure3, 
       filename = "figure3.png", dpi = 300, 
       type = "cairo", width = 25.5, height = 25.5, units = "cm")

# Table A1. Effects of Service Decline on Willingness to Coproducele A1 ####

stargazer(
  core.mod1,
  comp.mod1,
  copr.mod1,
  type="text", ci = TRUE, digits = 2,  star.cutoffs = c(0.05, 0.01, 0.001))


#Table A2. Interaction effect of service decline and core task encouragement####

stargazer(
  core.mod2,
  comp.mod2,
  copr.mod2,
  type="text", ci = TRUE, digits = 2,  star.cutoffs = c(0.05, 0.01, 0.001))

#Intext refs #####

#....effect of the service decline cue on willingness to engage in 
#overall coproduction (b = REF)
# B
fig2ci95$estimate[fig2ci95$dv == "Overall coproduction"] %>% round(digits=2)
# Cohens D  
fig2panbci95$Cohens_d[fig2panbci95$dv == "Overall coproduction"] %>% round(digits =2)

#.more apparent for core tasks, here we found an increased willingness of...
# B
fig2ci95$estimate[fig2ci95$dv == "Core tasks"] %>% round(digits=2)
# Cohens D  
fig2panbci95$Cohens_d[fig2panbci95$dv == "Core tasks"] %>% round(digits =2)

#... for complementary tasks
# B
fig2ci95$estimate[fig2ci95$dv == "Complementary tasks"] %>% round(digits=2)
# Cohens D  
fig2panbci95$Cohens_d[fig2panbci95$dv == "Complementary tasks"] %>% round(digits =2)

#....Z-test 

core.mod1
comp.mod1 


d<-coef(core.mod1)[2] - coef(comp.mod1)[2] # Difference in the betas

#coef(fit1)[2] # is the second coefficient in the model, i.e. the beta

var1<-summary(core.mod1)$coef[4]^2 

var2<-summary(comp.mod1)$coef[4]^2

ztest<- d / sqrt(var1+var2)

ztest %>% round(digits =2)

2*(1-pnorm(ztest)) %>% round(digits=3)


